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Abstract 


This  report  examines  the  problem  of  adaptively  tracking  a  maneuvering 
target  in  two  dimensional  space  using  passive  time  delay  measurements. 

The  target  is  free  to  maneuver  in  velocity  and  to  make  depth  changes  at 
times  unknown  to  the  observer. 

Tracking  is  accomplished  by  using  the  basic  polar  model  of  target  and 
observer  motion  previously  developed,  and  included  in  this  report  for 
the  convenience  of  the  reader.  However,  the  important  distinction  is  that 
now  a  nonlinear  prefilter  has  been  added  to  the  tracking  system.  This  leads 
to  two  major  benefits:  the  first,  is  that  the  need  for  extended  Kalman 
filters  is  completely  eliminated  which  gives  the  tracking  system  a  much 
larger  degree  of  ^robustness1^  than  it  previously  had.  The  second  benefit 
is  a  decoupling  of  the  depth  estimator  from  the  polar  range  estimator, 
which  considerably  reduces  the  computational  level  of  the  adaptive  tracking 
system.  ; 


TABLE  OF  CONTENTS 


CHAPTER 

I 

INTRODUCTION  AND  BASIC  TARGET  MODELING 

CHAPTER 

2 

TWO  DIMENSIONAL  SYSTEM  ANALYSIS  AND  DEVELOPMENT  OF  THE 
NONLINEAR  PREFILTER 

CHAPTER 

3 

STATISTICAL  ANALYSIS  OF  THE  LINEARIZED  RANGE  AND  DEPTH 
MEASUREMENT  ERRORS 

CHAPTER 

4 

DEVELOPMENT  AND  SIMULATION  RESULTS  OF  ADAPTIVE  DEPTH 
ESTIMATOR 

CHAPTER 

5 

DEVELOPMENT  AND  SIMULATION  RESULTS  OF  ADAPTIVE  RANGE 
ESTIMATOR 

CHAPTER  6 


CONCLUSION 


Chapter  1 


1.1  Introduction  to  Target  Tracking 

During  the  past  several  years  much  effort  has  been  spent  in  the 
development  of  sophisticated  digital  filtering  algorithms  for  tracking 
maneuvering  targets.  A  common  method  has  been  to  model  the  target 
dynamics  in  a  rectangular  coordinate  system  which  results  in  a  linear 
set  of  state  equations,  but  forces  the  measurements  to  be  nonlinear 
functions  of  the  state  variables.  With  this  model  an  extended  Kalman 
filtering  algorithm  is  frequently  used  both  to  provide  current  state 
variable  estimates  and,  by  a  one-step  prediction  process,  to  linearize 
the  next  measurement  vector.  This  method  works  moderately  well  until 
the  target  makes  an  abrupt  change  in  its  trajectory  in  response  to 
pilot  or  missile-guidance  program  commands.  In  this  situation  the 
velocity  and  position  estimates  can,  and  often  do,  diverge  from  the 
true  unknown  values.  The  inherent  problems  of  this  approach  can  lead 
to  large  bias  errors  and  sometimes  complete  filter  divergence. 

Earlier  work  on  the  maneuvering  target  tracking  problem  includes 
Jazwinski's  limited  memory  filtering  [1],  in  which  the  filter  gains 
are  prevented  from  decaying  to  zero.  Another  technique,  described 
bv  Thorp  [2],  involves  switching  between  two  Kalman  filters  in  res¬ 
ponse  to  a  detected  maneuver.  A  third  approach,  due  to  Singer  [3], 
models  the  target  trajectory  as  a  response  of  the  target  model  to 
a  time-correlated  random  acceleration.  With  this  method  additional 
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state  variables  are  used  to  generate  the  correlated  forcing  func¬ 
tions  which,  in  turn,  increase  the  dimension  of  the  Kalman  filtering 
algorithm.  In  this  manner  the  technique  provides  the  filter  with 
statistical  information  concerning  target  maneuvers  based  on  an  as¬ 
sumed  range  of  possible  aaelerat ions .  Singer's  method  was  subsequen¬ 
tly  extended  by  many  others. 

Parallel  to  the  effort  was  the  method  of  modeling  major  changes 
in  target  trajectories  by  a  semi-Markov  process.  An  application  of 
this  approach  to  tracking  maneuvering  targets  in  two-dimensions  by 
Moose  [4]  was  successfully  extended  by  Gholson  and  Moose  [5]  to 
three-dimensional  tracking. 

The  general  approach  which  uses  the  "adaptive  semi-Markov  man¬ 
euver  model"  of  [4]  and  [5]  implies  a  discretization  of  possible 
vehicle  accelerations  or  velocities.  The  estimation  algorithm  then 
views  the  maneuvering  vehicle  as  if  it  is  responding  to  commands  which 
are  modeled  by  a  semi-Markov  process,  i.'e.,  a  random  process  with  a 
finite  number  of  "states"  (commands)  which  are  selected  according  to 
the  transition  probabilities  of  a  Markov  process.  A  semi-Markov  pro¬ 
cess  differs  from  a  Markov  process  in  that  the  duration  of  time  in 
one  state  prior  to  switching  to  another  state  is  itself  a  random  var- 
ble  [6J.  Incorporating  the  semi-Markov  concept  into  a  Bayesian  estimator 
was  done  in  [4]  and  [5].  This  estimation  algorithm  provides  a  substan¬ 
tial  improvement  in  filter  stability,  which  means  that  large  bias 
errors  are  prevented  from  being  built  up  due  to  unmodeled  target  accel¬ 
erations.  An  important  aspect  of  this  adaptive  estimation  algorithm 
is  its  elimination  of  a  "growing  memory"  which  is  prevalent  in  many 


adaptive  filters. 


With  the  brief  history  of  the  maneuvering  target  tracking  pro¬ 
blem  presented  in  the  previous  section,  we  see  a  general  progression 
in  the  sophistication  of  tracking  filter  design  stemming  primarily 
from  the  method  in  which  the  unknown  target  accelerations  are  modeled. 
This  trend  is  graphically  outlined  in  Figure  1. 

Initially,  target  maneuvers  were  modeled  as  the  response  to 
uncorrelated,  zero-mean  variations  about  a  nonaccelerating  target, 
shown  in  Figure  1(A).  As  a  result,  the  estimation  algorithm  could 
follow  only  those  maneuvers  which  were  comparable  with  the  input  noise 
level.  Furthermore,  the  filtering  results  during  nonmaneuvering  sit¬ 
uations  were  seriously  degraded  due  to  the  uncorrelated  input  noise. 

As  shown  in  Figure  1(B),  Singer  [3]  attempted  to  model  large-scale 
maneuvers  by  assuming  a  time-correlated  input  process  and  incorpora¬ 
ting  the  statistics  into  the  subsequent  filter  design.  In  Figure  1(C) 
large-scale  target  maneuvers  were  modeled  as  a  stochastic  process  whose 
mean-value  switched  randomly  from  among  a  finite  set  of  predetermined 
values.  The  adaptive  estimation  algorithm  mentioned  in  the  previous 
section  could  then  be  used  to  track  the  maneuvering  target.  This 
method  was  seriously  restricted,  however,  by  the  requirement  of  a 
large  number  of  preselected  mean  values  in  order  to  ensure  convergence 
of  the  estimation  process.  It  has  been  shown  (6)  that  by  combining 
the  concepts  illustrated  in  Figure  1(B)  and  Figure  1(C)  the  number  of 
mean  values  required  to  prevent  filter  divergence  is  greatly  reduced. 
This  combination  is  illustrated  in  Figure  1(D).  The  primary  benefit 
of  this  approach  is  the  large  saving  in  computational  effort.  An 
additional  benefit,  at  least  from  a  subjective  viewpoint,  is  that  the 
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Figure  1.  Historical  development  of  maneuvering  target  model  inputs. 
(A)  Zero-mean  white  Gaussian  plant  disturbance.  (B)  Correlated,  zero- 
mean  plant  disturbance,  (C)  White  Gaussian  noise  with  randomly  switching 
mean.  (D)  Correlated  Gaussian  noise  with  randomly  switching  mean. 
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2.  Target  motion  model. 


time-correlated , randomly  switching, mean-forcing  function  more  ade¬ 
quately  models  real-world  target  maneuvers. 


The  basic  target  modeling  ideas  are  shown  in  Figure  2.  The  tar¬ 
get  trajectory  is  generated  by  the  random  selection  of  an  input  time- 
correlated  Gaussian  process  whose  mean  value  is  applied  to  the  tar¬ 
get  plant  dynamics  for  a  random  duration  of  time.  This  input  distur¬ 
bance  process  lasts  until  a  new  input  Uj  is  randomly  chosen  from  among 
a  finite  set  of  n  possible  inputs.  With  this  model  as  a  background 
and  using  an  appropriate  choice  of  state  variable  equations  to  repre¬ 
sent  target  dynamics,  either  submarine  or  aircraft,  it  is  possible 
to  develop  an  "optimal"  (in  the  minimum  mean-square  error  sense) 
tracking  filter  that  adaptively  learns,  then  quickly  adjusts  itself 
for  each  major  alteration  of  target  trajectory. 


1.3  Incorporation  of  Singer  Process  into  the  Target  Dynamics 
In  incorporating  the  correlated  process,  the  linearized 
polar  model  of  [8]  is  preserved.  To  this  end,  consider  a  target  whose 
motion  in  rectangular  coordinates  is  described  by 
x  =*  -ax  +  u  +  w' 


X  X 

•  I  _  t  , 

w  “  -  aw  +  w 


(1.3.1) 


X  XX 


where 


a  is  a  drag  coefficient 

ux  is  the  deterministic  input  in  the  x  direction  randomly 

chosen  from  a  set  of  N  possible  inputs. 

is  the  Singer  correlated  acceleration  process  acting 

in  the  x  direction  with  a  time  constant  i  «  1/a. 

c 


w^  is  a  white  Gaussian  random  process  acting  in  the  x  direction 


A  similar  set  of  equations  exists  for  the  y  and  z  directions. 
Defining 

■  x 


X2  ■  X 


X-,  *  w' 

3  x 


the  following  continuous  time  state  variable  model  is  obtained  for 
Equation  (1.3.1) 


0  10  xL  0 


0-al  x_+l  u  +  0 

2  x 


0  0  -a  x. 


(1.3.2) 


Discretizing  (1.3.2)  in  time  yields 
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0  E  F 
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x_  +  A  u  + 


w  (1.3.3) 
Xk 


where 


A  *  (1  -  e  aT)/a 

B  »  [1  +  (ae  aT  -  ae  aT)  /  (a  -  a) ]  /  (aa) 
C  -  (ctT  -  1  +  e‘aT)/a2 


D  *  [T  +  (aA  -  oJ)/(n  -  a)]/(aa) 


E  -  e 

F  *  (e  aT  -  e  aT) / (a  -  a) 


C  =  (J  -  A) /a  -  a) 
J  -  (1  -  e-aT) / a 


1.4  Development  of  the  Polar  Plant  Model 


(:<  “  O  ( y „  -  y„) 

:  s  o  ,•  •  ,  .  s  o  ,  •  •  , 

P  ■  -  (x  -  ::  )  f  -  (y  -  y  ) 

p  so  p  s  7c 


(1.4.2) 


p  ^  is  expanded  as  follows  keeping  only  the  linear  terms: 


3o  z  .  ,  So  ,  . 

-  -  )  +  - - ,  (v  -  y  ) 

s'k  "'k+l  Sk  ays :  k  Sk+1  Sk 


P,  ,,  *  0,  +  t+—  (x  -  x  )  +  r 
k+1  k  3x  ! ,  s  - 


+  4~ — |  (x  -  x  )  +  i  (y  -  y  ) 
3xolk  °k+l  °k  ayo!k  °k+l  °k 


(1.4.3) 


Assuming  a  linear  drag  model  for  the  Source  as  given  in  (1.3.1),  upon 
substituting  (1.4.1)  into  (1.4.3)  for  the  Source  connected  terms,  and  we  ge 
(xk+l  “  xk)  *  ^k  +  Bwxk  +  Cuxk 

(x  -  X  ) 

o  »  p  +  - j  [Ax  +  Bw'  +  Cu  +  Dw  ] 

k+1  k  p  'k  3  s^  s^  sx J  ;k 

(ys  "  V , 

+  -  ,  [Ay  +  Bw'  4-  Cu  +  Dv  ]  ■ 

o  ■  k  s  s  s  s 

y  y  y  'k 


(x.  '  Xn} 

s _ 0_ 

p  1  k  l"o 


Cx„  -  xi 


(y.  -  y  )  r„ 

3 _ O  I  l  < 


k+1  °k 


o,  '  o,  ' 

It  .<4-1  < 


(1.4.4) 

In  the  above  equation  and  in  ail  subsequent  analysis,  subscripts  s  and 
o  refer  to  the  Source  and  the  Observer,  respectively. 


Now 


and 
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x  ) 
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(1.4.5) 


(y. 


k+1 


o, 


y0.T 

1C 


Combining  terms  with  like  coefficients  in  (1.4.45  and  making  use  of 
(1-4.5), 
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(xs  '  X0) 


(ys  -  V 


ys  -  C 


(X.  ’  XJ 
.  s  o 

+  -  X 

P  O 


(ys  -  y0) 


?  O 

(1.4.8) 


(x  -  :.  )  (y  -  y  ) 

using  -  x  + - 

cop 


y  =  V  Cos  3 
1  O  0  so 


(1.4.9) 


where  VQ  =  observer  velocity,  3Sq  is  the  bearing  between  VQ  and  the  target. 


using  1.4.8  and  1.4.9  yields 

(xa  "  *  >  .  (y.  -  yJ  . 

— s - o_  +  — a - o_  v  .  $  *  v  Cos  6 

g  s  p  '  s  o  so 


(1.4.10) 


Substituting  Equations  (1.4.9;  and  (1.4.10)  for  the  coefficients  of 
T  and  A,  respectively,  in  (1.4.6)  and  collecting  like  terms  results  in 
the  following  state  variable  model: 

2  “  0  +  Ac,  +  3w 1  +  Cu  Dw  +  (A-T)  (V  Cost:  l1.  (1.4.11  ) 

k+1  k  s  s,  s„  o  so  < 

Pk  Pk  °k 


A  similar  approach  is  taken  to  develop  a  state  model  for  0^+^ 

By  assuming  chat  the  Observer  maintains  a  constant  velocity  for  long 

periods  of  time,  and  that  w1  is  a  zero  mean  correlated  GAUSSIAN  random 

sp 

process  acting  in  the  p  direction  the  following  state  equations  are 
easily  obtained. 


(1.4.12) 

The  underlying  constraint  for  the  state  model  (i.4.12)  requires  chat 
the  Observer  adopt  a  constant  velocity  profile. 


10 


Z  Channel 


Using  a  discretized  version  of  the  basic  linearized  drag  model  of 
Section  1.3,  and  letting  be  the  target  depth  the  following  state  model 
is  presented. 
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(1.4.13) 


The  state  model  is  also  subject  to  the  constraint  of  a  constant  velocit 
observer . 

It  is  interesting  to  observe  that  both  state  models  (1.4.12)  and 
(1.4.13)  are  identical  although  the  types  of  manipulations  involved  in  the 
derivations  are  quite  different.  For  example,  a  Taylor  series  expansion 
is  used  to  derive  (1.4.12)  whereas  no  such  expansion  is  used  in  deriving 


2.1  Two  Dimensional  System  Analysis 


To  develop  the  measurement  equations  use  is  made  of  Figure  (4)  below. 
Hassab  (7)  at  NUSC  developed  the  following  set  of  equations  for  the  dif¬ 
ference  in  arrival  times  and  Where  represents  the  time  delay  dif¬ 

ference  between  the  direct  path  along  (r)  and  the  surface  reflected  path. 


ocean  surface 


Figure  4.  Two  Dimensional  (Range  and  Depth)  Tracking  Geometry 
From  Hassab  we  have: 

Xl  «  [(r2  +  4  dQ2  -  4  dod)^  -  r]/C 

!  (2.1.1 

t2  =  [(r2  +  4  dk2  +  4  dkd)"*  -  r]/C 

2  2  2 

and  r  =  o  +  d 

C  »  speed  of  sound  in  water 
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larger  than  the  observer  depth  squared  (dQ“)  then  becomes,  where  p  Is 
the  polar  target  range. 


2  d  d 
o  T 

PC 


(2.1.2) 


If,  in  addition,  we  have  the  "shallow  water  ranging  situation,"  i.e.  that 
2  2 

(d,  / r  )  is  '<  1.  For  example,  if  r  _>  2.25  d  ,  our  error  of  approximation 


is  less  than  1.4 X,  thus  becomes 

2  d,  (d  -  dj 
k  w  T 


PC 


(2.1.3) 


where  d  =  total  ocean  depth, 
w 

Examining  the  set  of  equations  2.1.2  and  2.1.3,  we  see  we  have  two  non¬ 
linear  algebraic  equations  in  terms  of  our  state  variables  p  and  d^.  At  this 
point  in  the  past,  these  equations  were  linearized  in  a  Taylor’s  series  ex¬ 
pansion  to  yield  a  linear  measurement  of  p  and  d,j,.  By  doing  this  periodically 
in  time,  an  extended  Kalman  filter  tracking  system  (6)  was  developed  and 
reported  upon. 

It  was  decided  to  eliminate  the  extended  Kalman  filters  in  order  to 
reduce  system  complexity  and  computational  burden,  and  to  increase  the 
robustness  of  the  tracker  against  any  possible  filter  divergence  problems. 

This  was  done  by  developing  the  nonlinear  prefilter  discussed  in  the  fol¬ 
lowing  section. 


2.2  Development  of  the  Nonlinear  Prefilter 

In  order  to  linearize  the  time  delay  measurements,  we  notice  that  each 
is  a  nonlinear  function  of  the  system  state  variables  polar  range  and  target 
depth  (p,  d^) .  By  dividing  by  t and  letting  a^  =  (dQ/d^)  the  ratio  of 
observer  depth  to  depth  of  water  under  the  keel  we  have  the  following  expres¬ 
sion  containing  only  d_. 
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2.3  Measurement  Linearization 

Considering  first  the  case  of  target  depth,  define  the  random  process 
=  (zd^  -  d^,)  to  be  the  measurement  error.  Combining  equations  2.2.1, 

2.2.3  and  2.2.4  yields 
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(T1  +  V  dW 


T1d 
1  w 


(T1  +  aoT2)  +  (vl  +  aoV2}  (T1  +  aoT2} 


or 


a  d 

o  w 

tt 

\ 

V1X2  ~  V2T1 


TI  +  vl  +  aoV2 


(2.3.1) 


where  t  =  (t.  +  a  t„). 

T  1  o  2 

The  target  depth  error  term  is  thus  the  ratio  of  two  Gaussian  random 

processes  x  and  y.  The  terms  x  and  y  are  both  sums  and  differences  of 

the  zero  mean  Gaussian  processes  and  v Thus  they  are  both  strongly 

correlated,  and  in  the  case  of  the  denominator  y, nonzero  mean,  which  becomes 

very  important  in  determining  the  structure  of  the  density  function  of  v,. 

d 

From  Figure  (5)  we  see  the  structure  of  the  linearized  measurement  zd^  = 


(dT  +  V 


Figure  5.  Nonlinear  Prefilter  for  Target  Depth  Measurement 

We  now  want  to  examine  the  measurement  error  of  the  polar  range 
equation  =  p  +  .  Take  the  difference  of  equations  2.2.2  and  2.2.4, 
we  find  that 


15 


(2.3.2) 


(t  +  a  t,)  +  (v  +  a  v  )  (x.  +  a  x  ) 

1  O  L  1  o  L  1  o  l 


(vl  +  V2) 


vp  P  (i—  +  V,  +  a  v„) 


o  V 


where  x_  *  x,  +  a  x_  and 
1  l  o  / 


p  =  bQ/x^,,  the  true  polar  range. 


Equation  2.3.2  can  be  simplified  considerably  if  we  define  the  random  process 


to  be  (v^  +  '  ^us  vx  Gaussian,  with  zero  mean  and  variance  = 


2  2  2 

+  aQ  a2  '  ^he  measurement  error  as  shown  in  Figure  (7)  becomes 


-P  v_ 


v_  +  x„ 
T  T 


(2.3.3) 


Figure  6.  Nonlinear  Prefilter  for  the  Target  Range  Measurement 


1* 


P?(y)  is  given  by 


p  2  (y )  =  px[x  =  g-1(y)  1  dg-~^-v-- 


(3.1.1) 


Letting  x  =  g  (y)  *  -ky/ (v  +  a),  and  P^(x)  the  zero-mean  Gaussian  density 


results  in 


P0(y)  =  — — - 7  exp 

v/Tn  (y  +  a)“ 


o  2 

2o“(y  +  a) 


(3.1.2) 


By  making  the  substitutions  a  =  p,  k  =  t_,  x  =  v  ,  y  =  v^  the  measurement 

T  T  P 

2  2  2  2 

error,  and  o_  =  (o,  +  a  a  )  we  have  for  the  densitv  of  the  polar  range 

r  i  o  o 


measurement  error 


r  2  2-t 

-  t  v  ; 

T  o  { 

?  7 


P9(v  )  =  - -  exp  | - ~ - - 

oT/ 2tt  (p  +  v  )  [2J_  (p  +  v  )” 

T  p  u  T  p  - 


(3.1.3) 


The  nature  of  this  density  function  is  shown  in  the  following  figures  for 
the  target  and  observer  at  fixed  depths,  but  closing  in  range. 

It  appears  that  the  density  function  P7(v^)  given  by  equation  (3.1.3) 
is  of  a  form  that  does  not  appear  in  such  complete  sets  of  integrals  as. 

Tables  of  Integrals,  Series  and  Products,  by  GradshTevn  and  Rvzhik, 

Academic  Press,  1980.  In  addition,  it  does  not  seem  likely  that  there  is 

a  closed  form  for  the  expected  value  (v  =  -ax/'(x  +  k)),  or  for  the  expected 

P 

2  2  2  2  2 

value  of  =  a~xV(x  +  k)**  where  x  is  NCO^o4").  Thus  one  must  use  numerical 

P 

integration  in  order  to  get  an  approximation  for  the  mean  and  the  variance 
of  measurement  error  (Vp)«  Examining  the  density  function,  we  see  that 
it  is  a  smooth  continuous  curve.  As  target  range  closes  the  mean  and  variance 
becomes  smaller  and  the  density  becomes  more  and  more  symmetrical  approaching 
a  near  Gaussian  like  density. 

After  several  discussions  with  Dr.  A.  H.  Nuttall  of  NUSC,  New  London, 

CT,  it  was  decided  to  limit  any  measurements  of  c  which  exceeds  a  certain 
threshold.  This  would  guarantee  that  the  variance  of  the  measurement  error 


IS 


remain  bounded  and  constant  for  fixed  target  range.  This  is  due  to  the 

small,  but  finite  probability  (especially  at  long  ranges  of  p)  that  v^t^) 

equal  -  x^t^)  in  equation  (2.3.3)  at  some  time  t^.  In  other  words,  it  is 

possible  that  occasionally  a  bad  data  point  will  crop  up  and  must  be  dis- 

2 

carded.  It  was  found  that  only  a  few  stored  values  of  E(v  )  and  E(v  ) 

P  P 

were  needed  in  the  tracking  system  as  the  target  range  closed. 


3.2  Statistical  Analysis  of  the  Target  Depth  Measurement  Error 

Repeating  the  target  depth  measurement  error  equation  2.3.1,  we  have 


Vd  '  t 


a  d  ( v,  x0  -  v  x  ) 
o  w  12  2  1 


T  (tT  +  V1  +  V25 


where  v^  and  v^  are  the  statistically  independent,  zero  mean,  Gaussian  measure¬ 
ment  errors  associated  with  measuring  x^  and  X£* 

Letting  x  =  (av^  -  bv^)  and  y  =  (x^  +  v^  +  a^v^) ,  we  have  E(x)  =  0, 

2  22222  2  t  2  2 

E(x  )  -  °x  =  (a  flj  +  b  a,  ),  and  E(y)  =  xT,  cy  =  ?2  )  •  X°w 

define  (r)  to  be  the  correlation  coefficient  of  the  linearly  related  randon 
variables  (x)  and  (y) .  It  will  be  assumed  in  this  section  that  r  is  known 
and  will  be  developed  later  in  this  report.  The  joint  density  function 
P^(x,y)  for  correlated  Gaussian  random  variables  x  and  v  becomes 


P1(x,y)  =  a  exp 


-  6 


2rx (y  -  v)  (v  -  v)4-!"} 

•jj  2  ; 

XV  -> 

-v 


(3.2.1) 


where:  r  *  —  -X-^ 


x  -  y 


v  =  E(y) 

f,  c 

a  =  ^TT  ;x  ,y  x.  - 
3  =  (2(1  -  r2))"1 
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Now  let  z  equal  x/y  where  the  joint  density  function  P(x,y)  is  the  known 
Gaussian  distribution  of  (3.2.1).  In  order  to  find  the  density  function 
of  z,  we  need  to  identify  the  region  of  the  x,v  plane  that  is  valid.  To 


do  this,  let  z^  be  any  possible  value  of  z.  Then  for  the  probability  of 


z  “  |yj  —  z0’  we  have  two  regions  -  the  first  for  v  positive,  yields  x  <_ 
yz^,  and  the  second  for  v  negative  results  in  x  >  yz^.  Using  Figure  (8), 


the  probability  becomes  for  Region  I  and  II 


Region  II 


Figure  8.  The  Regions  of  Interest  for  Z  =  X/Y 


co  z  V  O  30 

PR[z  1  z o ]  =  /  /  °  P^(x,y)  dx  dy  +  /  /  P^x.v)  dx  dv  (3.2.2) 
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O' 

To  get  the  probability  density  functions  P(z^),  we  must  differentiate 

the  cumulative  probability  distribution  of  equation  3.2.2.  When  we  do 

this,  and  make  the  substitution  z  =  z ,  we  h3ve 

o 


P(z)  =  /  y  P ^ (vz ,  v)  dy  -  /  v  P  (yz.  v) 


(3.2.3) 
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where  P^(x  =  vz,  y)  is  the  joint  density  of  x,  v  with  vz  substituted  for  x. 
The  first  integral  ( I becomes 


=  a  j  v  exp  ' 


C, 


L 


2  2 

v  z  _  2ryz 


Jy 


(y  -  y)  + 


(y  -  v)2 


J  !  dv 


Expanding,  collecting  terms,  and  letting 


a  = 


8z  2rBz  8 


„ .  JSL  f„  .  I 

°x  3V  l  3 


(3.2.4) 


I  becomes:  ~ 
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Noting  that  /  ye 
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dt  and  that  the 


/2a 


right  hand  integral  is  the  standard  normalized  Gaussian  probability  integral 
b  ) 


1  -  F 
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we  have  for  I, 
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Likewise  for  I„  _2 
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(3.2.7) 


Combining  equations  3.2.3,  3.2.6  and  3.2.7,  we  have  for  the  density 
function  of  z 
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P(z)  =  f 
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where  a,  3  are  defined  in  3.2.1,  and  a,  b  in  equation  3.2.4. 

To  get  a  physical  feel  for  the  form  of  the  density  function  P(z)  where 

z  *  x/y,  if  we  let  E(y)  =  y  =  0  then  b  =  0,  and  P(z)  reduces  to  the  fol- 

2.-1 


lowing  when  g  is  replaced  by  (2  -  2r  ) 

Vy/TT? 


p(z)  «  f 
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r  x 
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+  ^x2(l  -  r2)1 


(3.2.9) 


r  x 


This  has  the  form  of  the  familiar  Cauchy  density  function  centered  at  z  = 

Jy 

Unfortunately,  in  our  case  y  =  ty  which  is  the  weighted  sum  of  time  delays 
t,  +  a  t„,  and  this  is  never  zero. 

1  O  2 

If  we  now  look  at  the  first  term  of  P(z)  in  3.2.8,  we  see  that  the 

_  2 

Cauchy  density  is  weighted  by  the  term  exp  j-  -^yj  .  Now  by  taking  the 
integral  of  the  entire  first  term  with  respect  to  z  we  have  just  the  previously 
mentioned  exponential  term  for  its  contribution  to  the  total  probability 
distribution  of  z.  If  this  contribution  is  large  (near  unity)  then  the 
second  part  of  equation  3.2.8  contributes  little,  and  p(z)  is  nearly  a 
Cauchy  density.  In  order  to  compute  this  weight  we  must  go  back  and  deter¬ 
mine  the  correlation  coefficient  r. 

3.3  Determination  of  Correlation  r  in  Target  Depth  Measurment  Error 

The  additive  measurement  error  v,  was  oreviously  found  to  be 

a 

a  d  (x0v  -  tv.) 
o  w  2  1  12 


TT  (tT  +  V1  +  aoV25 


When  we  originally  substituted  z  =  x/y,  we  let  x  be  a  Gaussian  process  given 

a  d 
o  w 

bv  -  (t.v,  -  t,v„),  and  y  also  Gaussian  is  given  by  fr  +  v  +  a  v  ).  We 

2  1  1  2  1102 

found  that  the  E(x)  =  0,  =  a“x^2  +  b2x  ^ t  E(y)  =  ty  and  z  ^  + 

7  7 

a  “a  “.  In  this  substitution  a  equals  a  d  and  b  equals  a  d 

o2  7  ow2T  owlT 

The  basic  definition  of  the  correlation  coefficient  r  is  given  by 
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r  * 


(3.3.1) 


txy)  -  t(x) 


where  the  expected  value  of  x  =  0.  The  first  numerator  term  becomes 


E(xy)  *  E[(av1  -  bv2) (tt  +  vL  +  a^)] 


2  2 
E(xy)  =  aa^  -  aQb 


(3.3.2) 


With  the  second  numerator  term  equal  to  zero  and  cx  and  cy  given  above 
the  correlation  coefficient  r  can  be  expressed  as 

f.  Vl  g22) 

*2  ®,2 

r  =  r  2  2-.1/2  2-,  1/2  (3.3.3) 

T1  ^  2  C2 
1  +  •—  1  +  a  ~ 

2,2  o  -  2 

l  t2  V  1  ' 


A  typical  scenario  might  be: 
p  =  25,000 


d  =  3,000 
w  ’ 

d  =  1,000 
o 

=  10  ms 

x2  =  80  ms 

a  =  1/2 
o 

then  if  a,  =  o  =>r  equals  .838 


if  a 2  =  2a2=>r  equals  .514 
indicating  a  strong  correlation  between  x  and  y. 


Chapter  4 


State  Estimation  and  Adaptive  Tracking  Svstem  Structure 


4.1  Introduction 

In  this  section,  we  will  discuss  the  basic  estimator  system  structure. 

We  will  make  use  of  the  linearized  measurements  containing  the  nonstationary, 
non-Gaussian,  statistics  that  was  previously  discussed  in  Chapter  3  of  the 
report . 

Referring  to  Figure  (9),  the  nonlinear  time  delay  measurements  Z 
and  are  fed  into  the  nonlinear  prefilter.  This  unit  develops  a  linearized 


Basic  Estimation  Structure 
Figure  (9) 


measurement  of  target  range  (p)  and  depth  (dt).  The  errors  in  measuring 
these  target  parameters  are  both  non-Gaussian  and  nonstationary  depending 
upon  the  geometry  of  the  tracking  situation.  As  target  range  closes,  or 
opens  the  mean  value,  and  variances  of  these  errors  change,  and  must  be  ac¬ 
counted  for. 

An  offline  computer  program  was  developed  to  generate  the  variation  in 

the  mean  value  and  variance  of  v  as  a  function  of  target  range  o.  Details 

P 

of  this  work  is  presented  in  the  accompanying  appendix  along  with  the  out¬ 
line  of  the  pitfalls  to  be  avoided  in  using  numerical  integration  routines 
on  functions  containing  an  isolated  singularlitv .  Results  shown  in  Figure 
(11)  illustrate  that,  for  stationary  Gaussian  measurement  errors  and  7, 

associated  with  time  delays  and  t9,  it  is  possible  to  simply  store  off- 

2 

line  as  a  small  catalog  the  necessarv  values  of  E(v  )  and  E(v  “) . 

P  P 

In  order  to  evaluate  the  mean  and  variance  of  the  target  depth  error 

v^  ,  data  was  generated  at  fixed  target  ranges.  It  was  then  operated  upon 
dT 

by  the  nonlinear  prefilter  shown  in  Figure  (10),  and  then  numerically 
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integrated  to  yield  a  set  of  means  and  variances  for  use  in  the  target  depth 
estimation  portion  of  Figure  (9) . 


4.2  Target  Depth  State  Estimation  and  Data  Generation 

A  conventional  Kalman  filter  was  developed  for  the  state  equation 

wzk  where  <!>  and  are  given  by  eq.  1.4.13  and  1.3.3.  The 

filter  of  the  form  =  4  x^  +  ^k+x  tZ^^Ck+l)  -  HS  x^]  gives  estimates  of 

the  target  depth  d^,(k+l)  which  is  the  first  or  upper  component  of  the 

estimated  state  vector  x,  , , . 

— k+1 

A  large  number  of  computer  runs  were  made  using  noisy  data  generated 

by  the  following  technique.  In  Figure  12  a  discrete  time  system  model  is 

presented  showing  the  development  of  noisy  time  delay  measurements  Z^  (k) 

and  (k) .  The  upper  portion  shows  the  generation  of  2^+x  =  ^  ^k  + ^ Uk 

+  1  W,  ]  where  <t,  I",  and  t 1>  are  discussed  in  eq .  1.4.12  with  cos  6  =1.0. 

The  deterministic  input  is  unknown  to  the  tracking  filter  and  serves  to 

generate  large  scale  target  maneuvers  in  velocity.  A  measure  of  randomness 

in  target  trajectory  (x^  =  o^  and  Xj  =  c^)  is  generated  by  applying  a  Gaussian 

random  input  to  the  simulated  target,  exponentially  correlated, forcing 

function  component  as  state  vector  x,  =  W  (k) . 

j  p 

Once  P^+x  generate<i.  it  is  acted  upon  in  a  nonlinear  manner  to 

generate  t,  and  t„,  which  when  added  with  the  Gaussian  random  measurement 

errors  v. ,  v_  produce  the  noisy  time  delav  measurements  Z  and  Z 

T1  T2 

With  the  background  completed,  we  will  now  present  some  computer  simula¬ 
tion.  The  results  shown  are  typical  of  those  that  were  observed  over  manv 
trial  tests.  In  Figure  13  (low  noise  case)  the  first  figure  is  a  plot  of 
the  linearized  nonGaussian  raw  data  Z^^(k).  The  target  is  at  a  fixed  range 


of  p  =  25,000,  and  at  a  mean  depth  of  600  feet  making  slow  random  changes 


o  =r 


2.5  msec  and 


on  the  order  of  +  30-40  feet.  Measurement  errors  of  * 
a,  *  5  msec  were  added  to  t  and  t~  in  Figure  12. 
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Noisy  Depth  Data  Z,  (k) 

a-p 

Figure  (13) 

In  Figure  14,  the  target  range  has  been  increased  to  40 ,000  and  o 
increased  by  a  factor  of  4.0.  Notice  the  data  shows  a  bias  due  to  an 
excessive  number  of  large  positive  readings.  In  the  future,  bad  data 
points  that  are  nonrealistic  will  be  "windowed  out." 


feet)  about  the  mean  value  of  600  whi 
Again  the  targets  range  is  fixed  at  2 
at  a  mean  depth  of  1000,  target  600, 


4.3  Depth  Estimation  for  Variable  Range  Targets 


The  work  presented  in  the  previous  section  dealt  only  with  targets  at  < 

i 

< 

fixed  ranges,  which  for  most  scenarios  is  a  very  unrealistic  situation. 

In  order  to  provide  depth  estimation  for  a  target  whose  range  is  time 
varying,  and  with  an  initial  estimate  of  range  and  depth,  a  study  was 
performed  to  analyze  the  effects  of  range  and  noise  variances  on  the  depth 
estimates  previously  developed. 

It  has  been  determined  that  both  range  and  depth  estimation  improve 
as  the  observer  increases  his  depth.  Now  with  target  depth  initially 
unknown,  and  this  being  an  exploratory  engineering  study  we  have 
decided  to  investigate  two  scenarios.  The  first  has  observer  at  1000  ft 
and  the  target  at  600,  the  ocean  depth  shallow  at  3000  ft.  The  second  is 
a  reversal  of  magnitudes  with  target  at  1000  and  observer  at  600.  Values 
obtained  then  would  allow  the  estimator  to  perform  anywhere  within  this 
range  without  going  to  a  new  set  of  tabulated  means  and  variances. 

What  we  need  as  the  target's  range  varies,  is  to  compensate  for  the 
changing  bias  introduced  In  the  nonlinear  data  operation  of  the  nonlinear 
prefilter  of  figure  (10)  page  26.  Now, the  dat3  generation  system  of 
figure  (12)  page  29  was  exercized  at  a  series  of  discrete  target  ranges 
=  5K,  10K,  15K,  ...,  90K.  Additive  Gaussian  noises  v^  and  vn  were 
generated  from  a  series  of  ten  independent  random  generators.  The  noisy 
measurements  zT^  and  zT7(o^)  were  then  fed  into  the  nonlinear  prefilter. 

Taking  the  output  and  subtracting  the  measurement  error  was  obtained. 

By  averaging  a  sequence  of  250  noisy  measurements  a  value  of  mean  and 
variance  of  v^  was  obtained,  for  each  target  range  c  .  This  was  repeated 
for  each  of  the  ten  random  sequences  to  produce  a  good  "Monti-Carlo"  set 


I 


of  mean  values  and  variances.  The  results  are  shown  in  figures  17  and  18. 
It  is  to  be  noted  that  the  curves  also  show  the  effect  of  increasing  the 
variance  of  the  Gaussian  measurement  errors  v^  and  V2  from  2  to  5  ms. 

Figure  19  gives  an  idea  of  the  form  of  the  density  function  of  the 
additive  noise  v^.  The  range  (p)  is  fixed  at  50K  and  the  scenario  of 
d  =  1000  and  depth  of  target  equal  to  600  is  used.  If  the  area  under 

OO 

the  curve  /  v,  P(v,)dv,  is  integrated,  a  value  very  close  to  the 
io  d  d  d 

desired  value  of  -50  feet  is  obtained. 
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Results  for  Depth  Estimation 

In  figure  20,  the  target  is  closing  the  range  from  p  =  80K  to 
20K  at  a  constant  depth  of  900.  An  initial  depth  estimate  5^(0)  equal  to 
200  was  chosen  with  the  observer  at  600,  and  ocean  depth  3000.  The 
standard  deviation  cf  3  ms  was  chosen  for  additive  noises  v^  and 
In  this,  and  the  figures  to  follow  raw  depth  measurements  zcjt  out  of  the 
nonlinear  prefilter  are  shown.  Note  the  decay  (on  the  average)  of  the 
noisy  measurements  as  the  target  closes  the  range,  making  the  received 
signal  to  noise  ratio  increase. 

Figure  21,  illustrates  the  case  of  higher  noise  power  =  a ^ 

=  5  msec  and  target  closing  from  50K  to  20K.  The  initial  depth  estimate 
is  again  chosen  to  be  200  with  target  at  an  unknown  fixed  depth  of  900. 

In  figure  22,  a  target  closes  from  10K  to  IK  and  makes  large-scale 
random  changes  in  depth.  The  estimate  is  initialized  at  dT(0)  =  200  and 
appears  to  follow  very  accurately.  Note  the  fact  that  the  unfiltered 
measurements  are  not  very  noisy.  This  is  due  to  the  close  range,  high 
signal  to  noise  ratio  that  is  in  effect  even  though  =  o  ^  =  3  msec. 

Figure  23  is  the  same  scenario  as  the  previous  figure,  but  with 
target  range  increased  from  50K  to  20K  and  depth  still  making  major 
random  variations.  Tracking  is  still  quite  good  with  the  exception  of 
the  tracker  lag  that  develops,  which  introduces  an  offset  of  100  oi  more 
feet.  This  ability  of  the  tracker  is  primarily  due  to  the  addition  of 
the  Singer  correlated  acceleration  which  was  previously  built  in  to  the 
filter.  It  is  also  to  be  pointed  out  that  the  tracker  continually  uses 
new  stored  values  of  the  mean  and  variance  of  the  nonGaussian  measure¬ 


ment  and  error  since  the  target  is  constantly  changing  range. 
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list  imute  vs  Time  for  Target  Closing  (80K  to  20K)  Range. 
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Figure  23:  Variable  Range  5Hk  to  ?0K  and  Variable  Depth  Target. 


A  final  sec  of  runs  were  made  to  study  some  different  scenarios 


that  might  be  of  interest  to  the  Navy.  The  first,  shown  in  figure  (24) 

is  a  tracking  situation  where  tKe  observer  is  maintaining  a  fixed  range 

of  p  »  30  K,  and  the  target  is  making  rapid  large  scale  random  depth 

changes.  The  filter  was  initialized  at  a  depth  estimate  of  200  feet 

for  the  target.  Target  depth  (dT)  is  shown  as  the  light  curve,  and  the 

depth  estimate  (d_)  as  the  heavy  curve.  The  raw  data  (z,_)  out  of  the 
i  a  r 

nonlinear  prefilter  is  also  shown.  The  worst  errors  were  on  the  order 
of  100  feet  lasting  for  about  30  time  samples,  or  300  sec  of  data. 

Figure  (25)  shows  the  target  at  a  fixed  range  of  p  =  30  K  proceeding 
at  a  constant  depth  of  250  feet.  The  observer  is  at  1000  maintaining 
the  30  K  fixed  range,  and  using  the  old  catalogue  set  of  means  and 

i 

variances  for  the  target  at  a  fixed  600  foot  depth.  A  fixed  bias  is 
observed  of  about  60  to  70  feet,  gradually  decaying  towards  zero  as 
a  new  set  of  means  and  variances  were  automatically  changed  in  the 
filter. 
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Chapter  5 

The  Polar  Range  Adaptive  State  Estimator 

The  heart  of  the  adaptive  filter  summarized  in  this  report,  is  in  the 
forming  of  the  total  estimate  of  the  target  states,  from  a  weighted  sum  of 

state  estimates  conditioned  on  the  N  possible  discrete  input  levels  u  ^ . 

P 

Consider  the  state  model  (1.4.12).  This  state  model  views  the  target  input 
acting  in  the  polar  direction  as  being  derived  from  a  time  correlated  Gaus¬ 
sian  density  having  a  mean  value  u^ .  Next  consider  a  series  of  N  such 
Gaussian  curves  with  displaced  mean  values  u  ,  i  =  1,  2,  ...,  N  and 
partially  overlapping  "tails"  as  shown  in  Figure  26.  If  a  bank  of  N  Kalman 
filters  is  formed,  each  filter  based  on  the  state  equations  of  equation 
(1.4.12)  with  the  deterministic  input  u^  being  a  different  one  of  these  N 


Series  of  N  Partially  Overlapping  Gaussian  Curves  with  Mean  Values  u 

Figure  (26) 

mean  values,  then  a  series  of  N  estimates  is  obtained,  each  conditioned  on 
a  different  Gaussian  curve.  Next  a  weighted  sum  of  these  estimates  is 


obtained  in  a  manner  to  be  disclosed  below,  and  this  weighted  sum  is  taken 
to  be  the  total  unconditioned  estimate  of  the  target  states  given  by  equation 
5.1. 

N  , 

x(k+l)  *  l  xU;(k+l)  W .  (k+1)  (5.1) 

i=l  1 

Now  as  the  target  executes  a  series  of  evasive  maneuvers  in  the  polar 
channel,  the  changing  input  necessary  to  produce  these  maneuvers  is  viewed 
as  randomly  switching  among  the  N  Gaussian  curves.  By  applying  semi-Markov 
statistics  to  this  switching  process  a  series  of  N  probabilities  VT,  i  = 

1,  2,  ...,  N  is  generated  where 


=  Pr  {target  input  is  being  derived  from  the  Gaussian 


curve  whose  mean  value  is  u 


(i) 


OR, 


W  (k+1)  =  Pr  (u  (k)  =  u 
i  0 


(i) 


Z(k+1) } 


(5.2) 


and 


x(l)(k+l)  =  E{x(k+1) ) u  (k)  =  u(l),  Z(k+1) } . 

o 

Equation  (5.1)  is  a  total  probability  expression  developed  from  the  basic 
relation  that 


x(k+l)  =  E{x(k+1) | Z(k+1) } 


is  the  optimal  mean-squared  estimate.  It  is  well  known  that  the  optimal 
input-conditioned  estimates  are  provided  by  suitably  matched  Kalman  filters. 
In  particular,  for  the  i—  filter 


-(i),,.^  _  (i)  fi  \  I  t.  ( i ) 

x  (k+1)  =  $x  (k)  +  fu 


+  KCk+I)  [z(k+l) 


where 


K(k+1)  -  M(k+1)  HT[HM(k+l)  HT  +  R] 
M(k+1)  =  $P(k)  <5T  +  ip  Q<pT 


and 


P (k+1)  =  [I  -  K(k+1)HJ  M(k+1)  . 


The  matrices  <f,  F  and  ip  are  used  to  denote  the  respective  coefficient 
matrices  in  (1.4.12). 

The  following  (repeated  here  for  completeness)  is  an  outline  of  the 
analysis  given  in  [5]  to  calculate  the  recursive  weighting  coefficients 
W^,  i  =  1,  2,  . ..,  N.  Defining  Z(k+1)  =  {Z(k),  z(k+l)},  apply  Bayes  Theorem 
to  (5.2)  and  we  obtain 


W.  (k+1) 


Pr  {u(k)  =  u^|Z(k)}  p{z(k+l)  !u(k)  =  ,  Z(k)'t 

p{z(k+l) 1 Z(k) } 


(5.3) 


The  time  varying  denominator  is  independent  of  i  and  is  therefore  common  to 
each  W^(k+1)  as  a  normalizing  constant.  The  first  numerator  factor  is  deter¬ 
mined  from  the  semi-Markov  input  process.  Expanding  this  factor  in  a  total 
probability  expression, 

N 

Pr{u(k)  =  u^^]Z(k)}  =  f  Pr{u(k)  =  u^^ju(k-l)  =  u^\  Z(k)}  W.(k) 

J=1  J 

And  since  Z(k)  has  no  influence  on  the  Markov  state  transitions. 


Pr{u(k)  =  u(l)|Z(k)}  =  7  8..  *..\(k) 

j=l  J1  J 


where  the  semi-Markov  probability  is 


(5.1) 


Ji 


=  Pr{u(k)  =  u(l);u(k-l)  =  u(j)> 


is  near  unity  for  j  =  i  and  near  zero  for  j  4  i. 


Combining  (5.3)  and  (5.4) 


W.(k+1)  =  C,  p{z(k+l)  |u(k)  =  t/L),  Z(k) }  £  9..  W.  (k) 

1  *  j=l  J1  3 

where  is  a  normalizing  constant 


(5.5) 


is  the  desired  recursive  relation  for  W^.  The  required  density  p  is  ap¬ 
proximately  normally  distributed  and  has  distribution 


P{z(k+1) |u(k)  =  u(l),  Z (k)  }  %  N{m.(k+1),  V(k+1)}  , 


(5.6) 


where 


and 


m.(k+l)  =  H[$  i(l)(k)  +  7  u(l)(k)] 


V (k+1)  =  [HM(k+l)  H  +  R] 


These  probabilities  can  thus  be  expressed  in  a  vector  recursive  form 
as 

Vi  ‘  S.  4  ,T  4  <5-;> 

where  is  a  time  varying  diagonal  matrix  whose  elements  have  been  previously 
computed  in  each  of  the  (N)  Kalman  filters  from  eq .  (5.6)  we  have 

P. .  =  N  [m.(k+l),  V(k+1)] 

The  term  C^.  is  a  normalizing  constant  computed  at  each  iteration  to 
ensure  that  the  sum  of  probabilities  equals  unity.  The  matrix  9  is  a 
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precomputed  matrix  whose  elements  contain  statistical  knowledge  of  the 
randomly  switching  plant  inputs.  In  practice,  the  diagonal  elements  are 
nearly  unity  and  the  off  diagonal  terms  are  set  equal. 

Although  it  might  appear  from  equation  (5.1)  that  an  entire  Kalman 
filter  algorithm  is  being  executed  N  times  (for  each  of  the  possible  inputs) 
at  each  time  iteration,  such  is  not  the  case,  since  the  process  and  measure¬ 
ment  covariances  Q  and  R  are  the  same  for  each  filter. 

What  differentiates  the  different  target  "states"  is  the  discrete 
levels  u^;  however,  the  target  dynamics  remain  unchanged!  The  entire 
covariance,  and  gain  analysis  of  the  Kalman  filter  algorithm  becomes  iden¬ 
tical  for  each  state  in  a  given  channel  and,  consequently ,  need  be  executed 
just  once  rather  than  N  times.  The  adaptive  filter  structure  is  shown  in 
Figure  27. 


5.2  Computation  of  Additional  Required  Covariance  Term 

After  an  initial  investigation,  it  has  been  found  necessary  to  com¬ 
pensate  for  the  lack  of  exact  statistical  knowledge  concerning  target 
maneuver  commands.  In  other  words,  the  u^  of  the  filter  bank  rarely, 
if  ever^  match  exactly  that  of  the  plant. 

The  mathematical  basis  for  the  error  compensation  technique  is 
based  upon  the  following  statistical  analysis. 

In  the  state  estimation  algorithm,  for  a  system  with  an  unknown  to 
the  tracking  filter  deterministic  input  u^,  and  random  disturbance  w^, 

x(k+l)  =  -*x(k)  +  Tu(k)  +  'fw(k)  (5.3) 

the  predicted  covariance  matrix  P(k+l/k)  is  given  by 

P(k+l/k)  =  t,P(k/kHT  +  fQ(k)?T.  (5.9) 
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However,  this  is  true  only  if  the  deterministic  input  to  the  plant  is 


known  exactly.  Consider  the  situation  when  a  mismatch  exists  between  the 
plant  input  u^  and  the  deterministic  input  u^  used  in  the  filtering 
algorithm.  Assuming  a  linear  estimator  of  the  Kalman  variety,  we  seek 
to  optimize  the  estimate 

x(k+l/k+l)  =  x(k+l/k)  +  K(k+l)[z(k+l)  -  Hx(k+l/k) ]  (5.10) 

where  x(k+l/k)  is  the  one  step  predicted  estimate  of  x.  The  gain  matrix 
K(k+1)  is  as  yet  unspecified. 

Let  x(k+l/k)  •  <f>x(k/k)  +  Fu* 

x(k+l/k+l)  =  [I  -  KH]  [<f>x(k/k)  +  Tu1]  +  KH]^  +  Tuk  +  Yw  ]  +  K  v  x 
where  K  =  K(k+1). 

Defining  the  estimation  error  at  time  (k+1)  as  =  -  x^+1  we  have 

x (k+1 /k+1)  =  [I  -  KH][*(x(lc/k)  -  x(k))  +  TCu1  ~  \)  -  *»k]  +  Kv(k+1)  (5.11) 

The  error  covariance  matrix  becomes 
P(k+l/k+l)  -  E(^+1  ^+1T) 

=  [I  -  KH]  [4>P(k/k)4>T  +  H'QH'T  +  D][I  -  KH]T  +  KRKT 

+  cross  terms  which  go  to  zero  under  the  expectation  operator.  The  new 
additional  quantity  D  has  the  value 

D  =  rE{(u1  -  uk)(ui  -  uk)T}rT  +  TElu1  -  uk)(i(k/k)  -  x(k))T}?T  (5.12) 

+  4>E {  (x(k/k)  -  x(k))(ui  -  uk):}rT. 

The  condition  for  choosing  K(k+1)  is  to  minimize 
trace  [P(k+l/k+l)] 

To  find  this  value  for  K(k+1)  it  is  necessary  to  take  the  partial 
derivative  of  P(k+l/k+l)  with  respect  to  K(k+1)  and  equate  it  to  zero. 

The  optimum  value  for  K(k+1)  is  found  to  be  of  the  standard  form 
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K(k+1)  =  P(k+l/k)  HT[HP(k+l/k)HT  +  R]_1  (5.13) 

where 

P(k+l/k)  =  [<j>P(k/k)<J>T  +  H'QY1  +  D]  (5.14) 

where  P(k+l/k)  is  computed  using  the  new  additional  error  term.  The 
subsequent  calculations  are  identical  with  those  for  the  Kalman  filter. 

Thus  the  effect  of  the  mismatch  in  the  two  inputs  is  propagated  throughout 


the  entire  Kalman  filter  algorithm.  From  earlier  work,  we  have  found  the 
estimates  of  the  target  states  to  be  close  enough  to  allow  the  assumption 

that  E{x(k/k)  -  x(k) }  -  0.  In  other  words,  our  estimater  is  an  unbiased 

i  i  T  T 

estimater.  Thus  D  reduces  to  TE{(u  -  u(k))(u  -  u(k)  )  }T  . 


Let  <5u  be  a  uniformly  distributed  random  variable  equal  to  (u^  -  u^) . 

Now  by  setting  the  range  of  the  random  variable  equal  to  -  Au  <  6u  <  +  Au 

2  2 

2 

the  variance  is  easily  shown  to  equal  Au  .  Here  Au  is  the  spacing  between 

12 

adjacent  pairs  of  levels  u^,  This  analytic  result  turns  out  to 

equal  very  closely  the  experimentally  obtained  values  in  the  1979  ONR 


Report  of  Moose  and  McCabe.  Thus  the  one  step  predicted  covariance  matrix 


P(k+l/k)  =  [d»P(k/k)d5T  +  w1  +  T  Au2  rT]. 

12 


(5.15) 
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In  order  to  have  a  benchmark  with  which  to  compare  the  adaptive  estimator's 
performance  the  following  scenario  was  devised.  A  target  at  the  extreme  range 
of  100,000  was  generated  on  a  closing  trajectory  of  25  ft/sec  at  a  depth  of  600. 

The  target's  plant  model  was  matched  exactly  with  that  of  the  filter.  This 
was  done  to  simulate  and  test  the  performance  of  the  optimal  (non-realizable) 
tracking  filter.  In  addition,  it  also  provided  a  test  of  the  nonlinear 
prefilter  and  the  adaptive  state  estimator  whose  weights  were  all  set  at  zero 
for  the  unmatched  filters  and  unity  for  the  single  matched  filter. 

Target  tracking  results  are  shown  in  figures  28  and  29  for  additive  noises 
of  1  and  5  m  sec,  respectively.  Examining  the  figures  We  see  that  for  each, 
two  different  simulations  are  shown.  They  represent  the  performance 
that  one  would  expect  out  of  the  unrealizable  filter  when  two  widely  different 
sequences  of  random  numbers  are  used  in  the  data  simulator  that  is  shown  in 
figure  12  on  page  29.  Notice  that  the  vertical  scale  has  been  changed  in 
figure  29  to  show  the  effect  of  the  nonlinear  prefilter.  The  range  estimates  (  p  ) 
are  shown  darker  than  the  actual  range  (  p  ) . 

The  data  simulation  of  Figure  12,  page  29  made  use  of  the  following 
parameters:  the  maneuver  time  constant  (1/a)  was  chosen  as  25  seconds,  the 
singer  correlated  acceleration  time  constant  (1/a)  in  equation  1.3.1  was 
selected  to  be  40  sec,  and  the  data  rate  was  chosen  to  be  one  sample  every 
10  seconds.  Although  studies  were  made  indicating  a  marginal  increase  in 
tracking  performance,  for  sample  intervals  of  2  and  5  seconds,  the  sample 
interval  was  retained  at  10  seconds  in  order  to  reduce  the  computation  burden 
for  real  world  applications. 

In  the  adaptive  estimator  of  Figure  (27),  6  levels  of  input  were 
chosen  to  span  the  expected  target  velocity  range  of  +30  (ft/sec)  for  an 
opening  target  and  -45  (ft/sec)  for  a  closing  target.  If  the  velocity 
ranges  were  greater  than  this  number  N, could  be  increased  slightly.  The 
u^'s  were  chosen  to  model  +30,  +15,  0,  -15,  -30,  -45  (ft/sec). 

The  next  figure  (30)  shows  the  realizable  adaptive  filter  tracking 
in  the  presence  of  high  noise  (5  m  sec).  The  target  makes  a  random  step 
change  In  velocity  at  time  k  =  150,  which  corresponds  to  1500  seconds  into 
the  sceanrio.  The  raw  data  out  of  the  nonlinear  prefilter  is  shown  almost 
filling  the  plot.  Notice  though,  that  the  variance  becomes  smaller  with 
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decreasing  range  then  grows  as  range  increases.  This  is  typical  of  all 
simulations  in  this  report,  and  point  out  the  nonstationary  random  processes 
involved.  Output  of  the  estimator  (darkline)  is  seen  to  fluctuate  about  the 
true  value  of  range  with  maximum  errors  on  the  order  of  5000  which  corresponds 
to  10-15%.  From  this  point  on,  in  the  presentation  of  the  results  of  the 
ranging  portion,  the  noisy  measurements  Zp  will  not  be  shown. 

In  figure  31,  we  see  the  results  of  tracking  a  medium  range  (0-45k) 
target  as  it  undergoes  a  maneuver  at  k  =  140  and  k  =  340.  The  target  is 
on  a  closing  trajectory  then  reverses  its  velocity  at  140  and  340.  The 
estimate  p  tracks  very  well  at  close  ranges  and  progressively  gives  noiser 
estimates  as  the  targets  range  increases.  This  is  due  to  the  variance  of 
the  measurements ,  which  are  not  shown,  increasing  with  range.  The  additive 
time  delay  noise  is  =  02  =  3  m  sec. 

Several  conclusions  can  be  drawn  from  figure  31.  The  first  is  that  we 
need  to  limit  ourselves  to  close  proximity  targets,  or  smooth  the  estimate 
0  (k),  or  operate  in  an  environment  with  smaller  additive  noises. 

Figure  32  shows  the  effects  of  decreasing  the  measurement  error  to  a 
standard  diviation  of  1  m  sec.  The  target  makes  a  maneuver  at  k  =  140  and 
the  estimate  (3  tracks  extremely  well  throughout  the  scenario.  The  next 
figure  (33)  shows  the  weighting  coefficients  ^  and  u.  as  they  vary  due  to 
the  target  maneuver.  It  is  seen  that  the  weights  switch  appropriately  from 
0  to  unity  as  required. 

The  last  figure  (34)  in  this  section  illustrates  velocity  estimation 
5  vs  time  for  a  target  that  abruptly  reverses  velocity  at  k  =  140.  The 
adaptive  state  estimator  very  quickly  maintains  track  as  shown.  The  additive 
time  delay  measurement  noise  has  been  increased  back  to  3  m  sec.  The  process 
noises  driving  the  plant  (target)  model  is  excessive  and  will  be  reduced  in  the 
following  work. 
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Target  Tracking  in  tlie  Event  of  Multiple  Maneuvers 
(Figure  31) 


SIGMA  EQUALS  1MSEC 


Target  Tracking  for  the  l.ow  Noise  Case 
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5-4  Modification  to  the  Adaptive  State  Estimator 


In  order  to  operate  in  a  noisy  environment  where  the  additive  measurement 
errors  attached  to  T;  and  are  greater  than  or  equal  to  3  m  sec,  it  becomes 
necessary  to  smooth  both  the  weighting  coefficients  and/or  the  output  estimate 
It  was  found  that  a  simple  first  order  digital  filter  of  the  form 

Vl  =  apk  +  bpk+l  (5'4-1) 

was  successful  for  good  smoothing  and  small  lag  in  tracking  maneuvers.  The 
term  represents  the  smoothed  output  range  estimate  at  t^+^  and 

the  "rough"  filter  input.  Coefficients  a  and  b  were  set  at  0.9  and  0.1 
respectively. 

In  addition,  to  output  smoothing  the  weighting  coefficients  w.  (k) 

i  =  1,  2,  6  were  averaged  as  follows: 

N-l 

Wf(k)  =  l  wt  (k-j)  (5.4.2) 

j  =  0 

Here  w^(k)  represents  the  time  averaged  output  of  the  moving  window  averager 
equation  5.4.2.  Numerous  simulations  were  performed  as  N  was  changed  from 
10  to  30.  The  results  follow. 

Figures  (35  and  36)  illustrate  the  filters  range  tracking  performance 
for  the  case  of  low  signal  to  noise  ratios  (long  ranges),  and  the  effect  of 
averaging  the  probability  weighting  coefficients  w^(k)  on  the  smoothness  of 
the  filters  output.  The  scenario  used  to  generate  figures  (35,36)  starts 
the  target  at  a  constant  range  of  50,000  feet  with  zero  closing  velocity  and 
then  changes  the  target  velocity  to  -15  ft/sec  after  2,500  seconds.  Figures 
35  and  36  show  the  filters  range  estimate  superimposed  on  the  target  range 
output  with  30  and  10  data  points  averaged  respectively.  It  is  easily  seen 
that  as  the  number  of  averaged  data  points  decreases  the  smoothness  of  the 
range  estimate  also  decreases. 

Figures  (37  and  38)  illustrate  the  performance  of  the  filter's  range  and 
velocity  estimates  and  the  adaptive  nature  of  the  weighting  functions.  The 
figures  also  show  how  the  number  of  averaged  data  points  w  ^(k)  effects  all  of 
the  above  mentioned  parameters.  The  scenario  used  to  generate  figures  (37,  38) 
starts  the  target  at  a  range  of  50,000  ft.  and  a  velocity  of  -30  ft/sec 
(-means  toward  observer)  and  changes  the  targets  velocity  to  +30  ft/sec  after 
2,500  seconds.  Figures  37  and  40  show  range  estimates  superimposed  over  the 
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State  Estimation  with  Output  Smoothing  and  (30)  Averaged  Weights 

( Figure  35) 


Estimate  of  Range  for  Approaching  /  Opening  Target  (20  points)  Averaged 

(Figure  37) 
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Estimate  of  Velocity  for  Maneuvering  Target  (20  points)  Averaged 

(Figure  38) 


plant  range  output  with  20  and  10  data  points  averaged  respectively.  Figures 
38  and  41  show  velocity  estimates  superimposed  on  the  plant  velocity  output 
with  20  and  10  data  points  averaged  respectively.  Figures  39  and  42  show  the 
learning  performance  of  the  weighting  factors  associated  with  the  filters  whose 
deterministic  inputs  (u^  most  nearly  matches  that  of  the  plot  output. 

Figures  39  and  42  average  20  and  10  data  points  respectively.  In  all  of  these 
types  of  curves  shown  in  figures  (37-42),  the  curves  become  smoother  as  the 
number  of  averaged  data  points  increases. 

Figure  43  shows  the  range  estimate  superimposed  over  the  target  range 
output  for  a  scenario  with  a  number  of  velocity  changes.  The  scenario  starts 
the  target  at  50,000  ft.  with  a  closing  velocity  of  -30  ft/sec.  At  1,500 
seconds  a  maneuver  is  made  with  velocity  changing  to  +30  ft./sec,  and  then 
again  changing  to  -30  ft/sec  at  3,500  seconds.  The  adaptive  tracking  filter 
quickly  determines  that  a  maneuver  has  been  made  and  adjusts  itself  to  provide 
a  good  estimate. 


Mu  I L i  Maneuver  Situation  with  Target  Undergoing  Maximum  Velocity  Changes 

(Figure  43) 


6  Conclusion 


An  adaptive  state  estimator  has  been  developed  and  extensively  tested  to 
track  a  target  making  random  large  scale  maneuvers  in  velocity  and  random  depth 
variations  as  well.  The  target/observer  scenario  is  constrainea  to  the  vertical 
plane  in  the  ocean  environment.  This  was  intentionally  done  so  as  not  to  compete 
with  well  established  bearing  tracking  programs. 

The  adaptive  estimator  made  use  of  a  nonlinear  prefilter  to  uncouple  the 
state  variables  that  model  target  motion  in  both  depth  and  range.  An  additional 
benefit  was  the  elimination  of  all  extended  Kalman  filters  in  the  tracking 
system.  This  results  in  a  more  vobust  tracker  and  significantly  fewer  computa¬ 
tions.  The  cost  of  doing  this,  is  that  the  linearized  measurements  contain 
nor.stationary  and  non-Gaussian  measurement  errors. 

System  inputs  to  the  tracking  system  consists  of  noisy  time  difference 
measurements  of  bottom/direct,  and  surface/direct  multipath  time  delays.  The 
adaptive  tracker  pre-filters  the  noisy  multipath  measurements  in  a  nonlinear 
operation  and  then  transmits  the  new  linearized  depth  and  range  measurements 
into  their  respective  filtering  channels.  The  depth  channel  gave  good  estimates 
as  the  target  underwent  random  depth  changes.  The  range  channel  was  more  complex, 
in  that  since  the  target  is  free  to  make  major  random  velocity  changes  it  required 
six  Kalman  filters  and  an  adaptive  weighting  technique  to  span  the  expected  range 
of  all  target  velocities.  Computationally  this  was  quite  easily  done  since  each 
filter  was  only  third  order,  and  had  the  same  Kalman  gain  and  covariance  matrix 
which  only  required  one  basic  computation  common  to  all  six  filters.  The  filters 
differed  only  in  the  deterministic  input  u^  i  =  i,  2,  ....  N  built  into  each. 

As  the  target  changed  range  in  5k  increments  new  means  and  variances  were  pro¬ 
gramed  into  the  filter  bank  thus  compensating  for  the  effect  of  nonstationarv 
linearized  measurement  errors. 

Overall  tracking  results  seem  quite  good,  especially  in  the  low  signal  to 
noise  ratio  cases.  An  on-going  effort  is  underway  to  study  techniques  of 
averaging  data  to  Gaussianize  the  measurement  errors  and  will  be  reported  upon 
in  a  later  report. 
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Appendix 


Problem:  Find  the  mean  and  variance  of  the  density  function  at  the  output 

2 

of  a  non-linear  filter  whose  input  is  N(0,  )  using  numerical  integration 

techniques. 


The  mean  and  variance  of  p(v^)  can  be  computed  using  the  following 


two  methods: 


A)  Calculate  the  output  density  function  from  the  input  which  produces 


p(v  )  and  integrate  as  usual. 
P 


2,2  mean  =  /  v  p(v  )  dv  = 
-T_  V  p  p  p 


p(v  ) 
p 


(v  +  p)^  a-  ^2tT  (v  +  p)2  .  r  ( 

d  r  T  d  variance  =  i  (v  -  m)  p( 


variance 


Note:  1)  The  tails  of  this  function  decrease  slowly  enough  to  cause  problems 


in  the  calculation  of  the  variance. 

2)  There  is  a  singularity  at  -p  which  must  be  dealt  with. 

These  problems  will  be  discussed  in  further  detail  latter  on. 

-p  vT 

B)  Let  —  ■  =  g(x)  and  calculate  the  mean  and  variance  as  follows: 

T  T 

CO 

A  "> 

mean  =  J  g(x)  p(x)  dx  where  p(x)  =  N(0, 


—  00 


Note;  Both  mean  and  variance  functions  have  singularities  at  -t^ . 


Calculation  of  Mean  and  Results 

The  mean  of  p(v^)  was  calculated  using  both  methods  A  and  B.  The 

numerical  integration  was  performed  by  the  IMSL  routine  DCADRE  which 

uses  cautious  Ronberg  extrapolation. 

The  singularity  in  eq.  A  caused  the  program  to  halt  due  to  exponent 

-99 

underflow  (i.e.  p(v  )  <  10  ).  To  avoid  this  problem,  p(v  )  was  defined 

P  P 

-99 

as  0  whenever  it's  value  was  less  than  10 

The  singularity  in  equ.  B  caused  the  program  to  halt  due  to  exponent 

99 

overflow  (i.e.  g(x)  p(x)  >10  ).  In  this  case,  it  was  assumed  that  the  area 

under  g(x)  p(x)  from  -(e  +  t„,)  to  iT  was  equal  to  that  from  to  -(t^  -  e) 

—6  * 

for  small  s  (e  <  10  ).  Since  these  two  areas  are  of  opposite  sign,  they 

should  cancel  each  other  allowing  integration  from  A  to  -(e  +  x^)  and  from 
-(xT  -  e)  to  B  without  affecting  the  final  answer.  (See  Fig.  4.) 

The  mean  was  calculated  using  dif¬ 
ferent  values  of  p  and  where  P*xT  » 


Fig.  4 


-  was  chosen  as  small  as  possible  without  causing  exponent  overflow. 


Calculation  of  the  Variance  and  Results 


The  variance  was  calculated  using  method  A. 


2  2 
v  7\ 


var 


V2?rT  (v  +  P)2  o  22 

e  dv 


•  (V  aT  /27 


F  (v) 


»  T 


p  Tt 


lim  r  *  T 

One  can  see  that  the  F(v)  =  -  e  =  const.  This  is  due  to 

V-KO  p- - 

Or*  </2u 


the  fact  that  lim  - — - —  =  ■—  =  1. 

V~  (V  +  P)2  V2 


So  at  large  values  of  |v|  the  integral  is  that  of  a  constant  which  becomes 
unbounded  as  the  limits  of  integration  increase  in  magnitude.  Therefore 
a  practical  integration  interval  of  -100  *P  to  100 *P  was  selected  to  produce 
a  bounded  and  useful  variance. 

The  variance  associated  with  each  P  are  listed  in  Table  1. 
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